function diff = delta_m_p(m,param)

    % derivative of delta(m) in the replacement region
    if param.zeta==0
        diff=param.sigma^2/2/(1-param.alpha).*1./(log(m)+param.B).^2.*(-1./m);
    else
        diff=(param.zeta+param.s.*param.lambda).*q_n(param.m_h,param.psi,param)./q_r(m,param.psi,param).^2.*(-q_r_p(m,param.psi,param));
    end

end